E95 brain analysis

library(Seurat)
library(ggplot2)
library(pheatmap)
library(tidyverse)

── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ tibble  3.1.3     ✓ stringr 1.4.0
✓ readr   2.0.1     ✓ forcats 0.5.1
✓ purrr   0.3.4     
── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

fig_dir = file.path("figs","brain_plot")
dir.create(fig_dir)

Warning in dir.create(fig_dir) : 'figs/brain_plot' already exists

srt95 = saveRDS(file = "data/srt_E95.Rds")

Error in saveRDS(file = "data/srt_E95.Rds") : 
  argument "object" is missing, with no default

srt33 = srt95[,srt95$orig.ident=="201104_33"]

library(dbscan)


db_out = dbscan(srt33@reductions$spatial@cell.embeddings,60)

ggplot(data=NULL,aes(x=srt33@reductions$spatial@cell.embeddings[,1],y=srt33@reductions$spatial@cell.embeddings[,2],col= factor(db_out$cluster)))+
         geom_point(size=0.8)+
  coord_fixed()+
  theme_classic()

sel_clu = names(table(db_out$cluster)[order(table(db_out$cluster),decreasing = T)][1])
sel_clu

[1] "1"

srt33 = srt33[,db_out$cluster==sel_clu]

saveRDS(srt33,file="data/srt33_clean.Rds")

sel_cluster = 23

DimPlot(srt33,reduction = "spatial",cells.highlight=colnames(srt33)[srt33$seurat_clusters==sel_cluster])+ggtitle(paste0("cluster ",sel_cluster))+coord_fixed()

srt_sel = srt33[,srt33$seurat_clusters==sel_cluster]

db_out = dbscan(srt_sel@reductions$spatial@cell.embeddings,60)

ggplot(data=NULL,aes(x=srt_sel@reductions$spatial@cell.embeddings[,1],y=srt_sel@reductions$spatial@cell.embeddings[,2],col= factor(db_out$cluster)))+
         geom_point(size=0.8)+
  coord_fixed()+
  theme_classic()

sel_clu = names(table(db_out$cluster)[order(table(db_out$cluster),decreasing = T)][1])
sel_clu

[1] "1"

srt_sel = srt_sel[,db_out$cluster==sel_clu]

library(edgeR)

Loading required package: limma

library(limma)

cnts = srt_sel@assays$RNA@counts
cnts = cnts[rowSums(cnts)>20,]
cnts = as.matrix(cnts)
allcounts = DGEList(counts=cnts)
allcounts = estimateDisp(allcounts, robust=TRUE)

Using classic mode.

  design_mat = model.matrix(~ srt_sel@reductions$spatial@cell.embeddings[,2])
  fit = glmQLFit(allcounts, design_mat)
  lrt = glmQLFTest(fit)
  top = topTags(lrt,n=Inf)
  top

Coefficient:  srt_sel@reductions$spatial@cell.embeddings[, 2] 

srt_markers = srt_sel@misc$markers
top10 <- srt_markers[srt_markers$p_val_adj<0.01,] #%>% group_by(cluster) %>% top_n(n = 30, wt = -p_val_adj) # %>%  top_n(n = 5, wt = avg_logFC)
top10$pct_diff = top10$pct.1-top10$pct.2
top10 = top10[top10$pct_diff>0.1,]
top10 = top10 %>% group_by(gene) %>% top_n(n=1,wt=avg_log2FC)

top10 = top10[top10$cluster %in% sel_cluster,]

top_sel = top$table[rownames(top$table) %in% top10$gene,]
top_sel = top_sel[top_sel$FDR<0.05,]
write.csv(top_sel,file =file.path(fig_dir, paste0("DE_markersforcluster_",sel_cluster,".csv")))
top_sel

sel_gene = c("Tcf7l2","Nr2f2","Wfdc1","Tal2","Boc","Zic4")
FeaturePlot(srt_sel,features = sel_gene,reduction = "spatial",coord.fixed = T,order=T,pt.size=0.5,ncol=3)
ggsave(file.path(fig_dir,paste0("markergene_",sel_cluster,"_region.pdf")))

Saving 7.29 x 4.5 in image

FeaturePlot(srt33,features = sel_gene,reduction = "spatial",coord.fixed = T,order=T,ncol=3,pt.size=0.01)

ggsave(file.path(fig_dir,paste0("markergene_",sel_cluster,"_allpuck.pdf")))

Saving 9 x 4 in image

sel_clusters = c(23, 27)
srt_sel = srt33[,srt33$seurat_clusters %in% sel_clusters]

db_out = dbscan(srt_sel@reductions$spatial@cell.embeddings,60)

ggplot(data=NULL,aes(x=srt_sel@reductions$spatial@cell.embeddings[,1],y=srt_sel@reductions$spatial@cell.embeddings[,2],col= factor(db_out$cluster)))+
         geom_point(size=0.8)+
  coord_fixed()+
  theme_classic()

sel_clu = names(table(db_out$cluster)[order(table(db_out$cluster),decreasing = T)][1])

srt_sel = srt_sel[,db_out$cluster==sel_clu]


cnts = srt_sel@assays$RNA@counts
cnts = cnts[rowSums(cnts)>20,]
cnts = as.matrix(cnts)
allcounts = DGEList(counts=cnts)
allcounts = estimateDisp(allcounts, robust=TRUE)

Using classic mode.

  design_mat = model.matrix(~ srt_sel@reductions$spatial@cell.embeddings[,2])
  fit = glmQLFit(allcounts, design_mat)
  lrt = glmQLFTest(fit)
  top = topTags(lrt,n=Inf)


srt_markers = srt_sel@misc$markers
top10 <- srt_markers[srt_markers$p_val_adj<0.01,] #%>% group_by(cluster) %>% top_n(n = 30, wt = -p_val_adj) # %>%  top_n(n = 5, wt = avg_logFC)
top10$pct_diff = top10$pct.1-top10$pct.2
top10 = top10[top10$pct_diff>0.2,]
#top10 = top10 %>% group_by(gene) %>% top_n(n=1,wt=avg_log2FC)

top10 = top10[top10$cluster %in% sel_clusters,]

top_sel = top$table[rownames(top$table) %in% top10$gene,]
top_sel = top_sel[top_sel$FDR<0.05,]
write.csv(top_sel,file =file.path(fig_dir, paste0("DE_markersforcluster_23_27.csv")))
top_sel

sel_gene = c("Tcf7l2","Dmrta2","Rprm","Tubb3","Pcdh8","Shh")
FeaturePlot(srt_sel,features = sel_gene,reduction = "spatial",coord.fixed = T,order=T,pt.size=0.5,ncol=3)
ggsave(file.path(fig_dir,paste0("markergene_23_27_region.pdf")))

Saving 7.29 x 4.5 in image

FeaturePlot(srt33,features = sel_gene,reduction = "spatial",coord.fixed = T,order=T,ncol=3,pt.size=0.01)

ggsave(file.path(fig_dir,paste0("markergene_23_27_allpuck.pdf")))

Saving 9 x 4 in image

srt_sel$MHB = "no"
srt_sel$MHB[srt_sel$ycoord>2900 & srt_sel$ycoord < 3180 & srt_sel$xcoord>4900] = "yes"

DimPlot(srt_sel,reduction = "spatial",group.by = "MHB")+coord_fixed()

sel_gene = rownames(markers.MHB[markers.MHB$p_val_adj<0.05,])
FeaturePlot(srt_sel,features = sel_gene,reduction = "spatial",coord.fixed = T,order=T,pt.size=0.3,ncol=4)

ggsave(file.path(fig_dir,paste0("markergene_23_27_MHB.pdf")))

Saving 10 x 8 in image

eye genes

eye_genes = c("Rax", "Vax1", "Vax2", "Pax2", "Six6")
FeaturePlot(srt33,features = eye_genes,reduction = "spatial",coord.fixed = T,order=T,ncol=3,pt.size=0.01)

srt33_eye = ScaleData(srt33,features = unique(c(VariableFeatures(srt33), eye_genes)) ,verbose=F)

expand_genes = c()
for(ge in c("Rax","Vax1","Six6")){
  print(ge)
  corr = t(cor((srt33_eye@assays$RNA@scale.data[ge,]), t(srt33_eye@assays$RNA@scale.data )))[,1]
  corr = corr[!is.na(corr)]
  corr = corr[order(corr,decreasing = T)]
  expand_genes = c(expand_genes, names(head(corr))[2:6])
}

[1] "Rax"

Warning in cor((srt33_eye@assays$RNA@scale.data[ge, ]), t(srt33_eye@assays$RNA@scale.data)) :
  the standard deviation is zero

[1] "Vax1"

Warning in cor((srt33_eye@assays$RNA@scale.data[ge, ]), t(srt33_eye@assays$RNA@scale.data)) :
  the standard deviation is zero

[1] "Six6"

Warning in cor((srt33_eye@assays$RNA@scale.data[ge, ]), t(srt33_eye@assays$RNA@scale.data)) :
  the standard deviation is zero

srt33_eye = FindNeighbors(srt33_eye,dims = 1:5,verbose=F)

sel_spa = srt33_eye@reductions$spatial@cell.embeddings[srt33_eye$seurat_clusters==1,]
plot(sel_spa)

db_out = dbscan(sel_spa,eps=50)

ggplot(data=NULL,aes(x=sel_spa[,1],y=sel_spa[,2],col=as.factor(db_out$cluster)))+
  geom_point()

idx = rownames(sel_spa)[db_out$cluster==1]
srt33_eye$eye = "no"
srt33_eye$eye[colnames(srt33_eye) %in% idx] = "yes"

paste(sum(srt33_eye$eye=="yes"),"/",ncol(srt33_eye[,srt33_eye$RNA_snn_res.1 %in% c(0,1,15,21,23,27)]),"=",sum(srt33_eye$eye=="yes")/ncol(srt33_eye[,srt33_eye$RNA_snn_res.1 %in% c(0,1,15,21,23,27)]))

[1] "105 / 4824 = 0.0217661691542289"

paste(sum(srt33_eye$eye=="yes"),"/",sum(srt33_eye$RNA_snn_res.1==0),"=", sum(srt33_eye$eye=="yes")/sum(srt33_eye$RNA_snn_res.1==0))

[1] "105 / 1499 = 0.0700466977985324"

paste(sum(srt33_eye$RNA_snn_res.1==0),"/",ncol(srt33_eye),"=", sum(srt33_eye$RNA_snn_res.1==0)/ncol(srt33_eye))

[1] "1499 / 18260 = 0.082092004381161"

save.csv(marker.eye[marker.eye$p_val_adj<0.01 & marker.eye$avg_log2FC>0,],file="data/eyes_marker.csv")

Error in save.csv(marker.eye[marker.eye$p_val_adj < 0.01 & marker.eye$avg_log2FC >  : 
  could not find function "save.csv"

expand_genes_eye = c()
for(ge in c("Cp","Vwc2")){
  print(ge)
  corr = t(cor((srt33_eye_sel@assays$RNA@scale.data[ge,]), t(srt33_eye_sel@assays$RNA@scale.data )))[,1]
  corr = corr[!is.na(corr)]
  corr = corr[order(corr,decreasing = T)]
  expand_genes_eye = c(expand_genes_eye, names(head(corr,n=10))[2:10])
}

[1] "Cp"
[1] "Vwc2"

DimPlot(srt33_eye_sel,reduction = "spatial",pt.size = 4)+coord_fixed()
ggsave(file.path(fig_dir,"eye_subset_clustering.pdf"))

Saving 7.29 x 4.5 in image

DotPlot(srt33_eye_sel,features = unique(c(unique(top10$gene),"Cp","Vwc2")),cols="Spectral")+coord_flip()
ggsave(file.path(fig_dir,"eye_subset_markers_dotplot.pdf"))

Saving 7.29 x 4.5 in image

FeaturePlot(srt33_eye_sel,features = c("Vwc2","Jkamp","Cp","Tlk1","Pon2","Ict1"),coord.fixed = T,sort.cell = T,reduction = "spatial",ncol = 3)

Warning: The sort.cell parameter is being deprecated. Please use the order parameter instead for equivalent functionality.

ggsave(file.path(fig_dir,"eye_subset_markers_spatial.pdf"))

Saving 7.29 x 4.5 in image

pp = FeaturePlot(srt33_eye,features = rownames(marker.eye)[1:12],order = T,raster = F,ncol = 4,coord.fixed = T,reduction = "spatial",max.cutoff = "q99")
ggsave(file.path(fig_dir,paste0("spatial_eye_markers.pdf")),plot = pp,width = 14,height = 9)

srt36 = srt95[,srt95$orig.ident=="201104_36"]

pp = FeaturePlot(srt36,features = rownames(marker.eye)[1:12],order = T,raster = F,ncol = 4,coord.fixed = T,reduction = "spatial",max.cutoff = "q99")
ggsave(file.path(fig_dir,paste0("spatial_eye_markers_36.pdf")),plot = pp,width = 14,height = 9)

load("~/scRNAseq_reference_datasets/WT.Robj")

FeaturePlot(WT,features = c("Six6","Cp","Vwc2"),order=T,split.by = "stage")

Warning in FeaturePlot(WT, features = c("Six6", "Cp", "Vwc2"), order = T,  :
  All cells have the same value (0) of Six6.
Warning in FeaturePlot(WT, features = c("Six6", "Cp", "Vwc2"), order = T,  :
  All cells have the same value (0) of Cp.
Warning in FeaturePlot(WT, features = c("Six6", "Cp", "Vwc2"), order = T,  :
  All cells have the same value (0) of Vwc2.

WT$eye_cluster = WT@assays$RNA@data["Six6",]>0.7
table(WT$eye_cluster)

FALSE  TRUE 
88322   457 

sparse.cor2

table(WT$eye_cluster,WT$stage)

       
        WT_65 WT_70 WT_75 WT_80 WT_85
  FALSE  2280 13218 27417 23298 22109
  TRUE      0     0     5    46   406

table(WT$eye_cluster,WT$cluster)

       
           0    1   10   11   12   13   14   15   16   17   18   19    2   20   21   22   23   24   25   26   27   28   29    3   30   31   32   33   34   35   36   37   38   39    4   40
  FALSE 3860 4211 3638 2578 2753 2077  862 1907  956  787 2255 5096 2908 2386  887 3041  298 3708 1064 2857  290  530 2650 2322 1819 1594 1158 1297 1555 1972 1606 3137  338 2575 1081 1434
  TRUE     0   85  168    1    2    0    0    0    1    0    0    0    0    0    0    0    0  175    0    1    0    0    0    0    0    1    0    0    2   16    1    0    1    1    0    0
       
          41    5    6    7    8    9
  FALSE  482 3635 2509 1311 4808 2090
  TRUE     0    0    1    0    0    1

WT$eye_clu_fine = "no"
WT$eye_clu_fine[ WT$cluster==10]="clu10"
WT$eye_clu_fine[WT$eye_cluster & WT$cluster==10]="eye_10"
WT$eye_clu_fine[ WT$cluster==24]="clu24"
WT$eye_clu_fine[WT$eye_cluster & WT$cluster==24]="eye_24"
table(WT$eye_clu_fine)

 clu10  clu24 eye_10 eye_24     no 
  3638   3708    168    175  81090 

markers_24 = FindMarkers(WT,ident.1 = "24",group.by = "cluster",only.pos = T,min.diff.pct = 0.1,verbose = F)

eye.24_ref = FindMarkers(WT,ident.1 = "eye_24",ident.2 = "clu24",group.by = "eye_clu_fine",min.diff.pct = 0.1,verbose = F)

eye.24_ref = eye.24_ref[eye.24_ref$p_val_adj<0.01,]
eye.24_ref = eye.24_ref[eye.24_ref$pct.1>0.3 & eye.24_ref$pct.2<0.1,]
eye.24_ref

srt85 = readRDS(file = "data/srt_E85.Rds")

FeaturePlot(srt85,features = c("Six6","Aldh1a3","Foxg1","Cp","Vwc2"),split.by = "orig.ident",ncol=4,reduction = "spatial",coord.fixed = T,order=T)

Warning in FeaturePlot(srt85, features = c("Six6", "Aldh1a3", "Foxg1", "Cp",  :
  All cells have the same value (0) of Six6.
Warning in FeaturePlot(srt85, features = c("Six6", "Aldh1a3", "Foxg1", "Cp",  :
  All cells have the same value (0) of Six6.
Warning in FeaturePlot(srt85, features = c("Six6", "Aldh1a3", "Foxg1", "Cp",  :
  All cells have the same value (0) of Aldh1a3.

ggsave(file.path(fig_dir,"eye_markers_E85.pdf"))

Saving 46 x 10 in image